Spatially Dependent Parameter Estimation and Nonlinear Data 



Given multiple images that describe chaotic reaction-diffusion dynamics, parameters of a 
PDE model are estimated using autosynchronization, where parameters are controlled by syn- 
chronization of the model to the observed data. A two- component system of predator-prey 
reaction-diffusion PDEs is used with spatially dependent parameters to benchmark the meth- 
ods described. Applications to modelling the ecological habitat of marine plankton blooms by 
nonlinear data assimilation through remote sensing is discussed. 

1 Introduction 

Parameter estimation in ODEs and PDEs has developed into a vast field in applied mathematics 
and control engineering. For models representing important physical processes, accurate estimates 
of appropriate model parameters may help inform short-term management decisions by model 
forecasting. However, to forecast a system one requires not only accurate parameter estimates, 
but also full knowledge of the initial state of the system. There are widely varying and powerful 
methods for parameter estimation of spatio-temporal systems including, but certainly not limited 
to Kalman filter methods [13, 1, 18], multiple shooting methods [7, 6], and adjoint methods [8]. 
A method of parameter estimation based on synchronization has drawn substantial interest [11, 
9, 10, 17, 20, 19, 16, 12, 15, 2, 14]. Applications include communications and cryptography [17], 
electronics and circuit dynamics [10, 4], and cardiac cell dynamics [15] to name just a few. There 
are currently several methods to estimate parameters based on synchronization. One approach is 
to optimize a time-averaged synchronization error on which synchronization acts as a regularizing 
force; the optimization problem of finding the minimum synchronization error in parameter space 
is well-posed [10, 12, 15, 2, 14]. Our interest here will be based on an approach to force a response 
model to adapt to observed data by developing additional equations for the parameters that depend 
on the synchronization error [9, 16]. 

To estimate model parameters by synchronization, we exploit a special variation of synchroniza- 
tion called "autosynchronization" . For systems of ODEs, an observed scalar time series is coupled 
to a response system during model simulation. The goal of this feedback is to cause the response 
system to synchronize to the drive system. Ideally a proof of convergence follows by demonstration 
of an appropriate Lyapunov function [9, 20]. In [19, 16] we see some generalizations of how to derive 
synchronization schemes for many systems including the case where, apriori, we do not know the 
model form of the drive system [16]. By autosynchronization, we can recover the model parameters, 
the current model state, and in some cases, a model form for an observed system. 

Stating an autosynchronization problem in the ODE setting, we require a drive system 
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u t = f(u,p), 



(1) 



from which we are able to sample data with (unknown to us) parameters p G R m . Then we must 
state a response system 

v t =g(u,v,q) (2) 

with the same model form as the drive system if q = p. By "same" we mean in as far as possible 
by our understanding of the underlying physics. Then the goal is that when u is coupled forward 
into Eq (2), then Eq (2) will synchronize with Eq (1) and u — > v. Furthermore, parameter ODEs 
are given by 

q t =h(u,v,q) (3) 

so that (v, q) — > (u, p) as t — > oo. 

The idea of synchronization was extended to one-dimensional systems of PDEs in [4] and two- 
dimensional systems in [2], where the authors considered the Grey-Scott and Barkely reaction- 
diffusion systems respectively. In these works, the authors observed synchronization of an infinite- 
dimensional system by coupling the drive and response systems at only a finite number of grid 
points. Further work examines parameter estimation for given PDE systems using optimization over 
the synchronization error surface as discussed above. The authors observe single-species assimilation 
as they drive the PDE system to synchronization while coupling with only one species [2]. However, 
none of these utilize autosynchronization. 

In many systems, it is very reasonable to expect that model parameters need not be spatially 
homogeneous. For example, taking our problem of interest, spatial inhomogeneity in parameter val- 
ues may be central when constructing models for coastal algal blooms, since plankton growthrate 
is affected by near-shore nutrient runoff and upwelling [5]. More to that point, ocean fronts and 
eddies cause flow-induced long-term inhomogeneities in the ocean which results in a formidable 
spatial structure for density profiles in the ocean [5]. Whether inhomogeneities be the result of 
the flow dynamics or of "boundary conditions" from nutrient runoff, they are an important con- 
sideration for modelling ecology over large coastal domains. Thus it is reasonable to argue that a 
biophysics-based model should accept spatially dependent parameters. 

A drawback to the aforementioned methods, including Kalman filter and multiple shooting 
methods, is that they do not consider spatially dependent parameter values, a priority noted in [13]. 
Parameter estimation by filtering methods adapted for PDEs can be computationally expensive [6] . 
Furthermore, we have found that some filtering methods have trouble during periods of exponen- 
tial growth, such as might be expected during plankton blooms. Optimizing the time-averaged 
synchronization error in some function space is far more complicated than the finite-dimensional 
alternative with scalar parameters as in [2]; optimization methods may not be practical. 

Our work aims to extend the method of parameter estimation for PDE systems by synchro- 
nization to autosynchronization, especially including autosynchronization with spatially dependent 
parameters. Thus, we investigate observed data from the PDE drive system 

u t (x,y,t) ={(u(x,y),p(x,y)) (4) 

with parameters p(x,y) 6 C°(f2) and a response system 

v t (x, y, t) = g(u(x, y),v(x, y), q(x, y)). (5) 

We formulate an associated system of PDEs for the parameters of Eq (5) 

1t(x,y,t) = h(u(x,y), v(x,y)) (6) 

with the goal that (v,q) — > (u, p) as t — > oo. We design our methods considering a benchmark 
system of reaction-diffusion PDEs. Since we know the model form of the drive system and the 
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parameters used to build the observed data, we can compare our estimated parameters with the 
exact parameters. 

In preview of the paper layout, we begin by introducing the reaction-diffusion equations that we 
will use as the drive system. We discuss how this system is solved numerically and the parameters 
used to simulate complex spatiotemporal dynamics. Next, we implement the response system 
and show the parameter PDEs used to find autosynchronization. We demonstrate the power of 
the systems of PDEs to autosynchronize by employing three different spatial functions for the 
parameters. Next, we show the estimated parameters and the convergence plots for both state 
variables and parameters to the correct values. Finally, we give an improvement on the response 
system that admits autosynchronization wherein only one species is sampled, which is an important 
breakthrough for applications since generally only the phytoplankton is easily observable. 



2 The Parameter Estimation Method 

Consider the system of two PDEs as given in [5] , 



BP P7 

_ = AP + P(1-P,-— , (7) 
dZ PZ 

m = AZ + k ¥Th- mZ > 

on a compact connected two-dimensional domain, J7, with zero- flux boundary conditions. 

In terms of the biology of the model, the system represents a dimensionless reaction-diffusion 
model for phytoplankton-zooplankton predator prey dynamics in a horizontal layer where vertical 
distributions of plankton are considered uniform. For simulations shown here, we choose to be 
a rectangle of size 864 x 288. Although shown in dimensionless form, the model is derived from 
principles in which phytoplankton concentrations obey a logistic growth and are grazed upon by 
zooplankton, following a Holling-type II functional response. First classified by Holling, [3], the 
Holling-type II functional response assumes a decelerating growth rate such that the predator, 
or consumer, is limited by its ability to efficiently process food. Zooplankton g row at a rate, k, 
proportional to phytoplankton mortality and are subject to a natural mortality rate m. For certain 
parameters, this system gives rise to transient spiral pattern behaviour on its way to spatially 
irregular patchy patterns [5]. In [5], parameters are set to k = 2, h = 0.4, and m = 0.6. For 
homogeneous initial plankton distributions, the system remains in a homogeneous state for all time 
so we use the perturbed initial conditions found in [5]. 

We solve this system using a finite difference method, with a nine-point center difference stencil 
for spatial derivatives and forward Euler time stepping. Our simulations use spatial discretization 
with dx = 2 and and Euler time step of dt = .2. The model output is treated as an image sequence 
given by a particular (known) model form but with unknown parameters k and m, to be determined. 
Thus we will mimic our target application of remote sensing oceanographic images of hyperspectral 
images filtered to reveal plankton blooms. Further, we will be interested to allow k and m to vary 
spatially as functions, k(x,y) and m(x,y). 

Our interest in this PDE model stems from our work in remote sensing, to build a better un- 
derstanding of our ocean's ecology. Particularly, we aim to predict short term behavior of coastal 
algal blooms. Such a system may in principle be modelled by estimating parameters directly from 
observed data in the field. However, hyperspectral satellite imagery provides the observed data 
to which we would synchronize a response model in hopes of autosynchronization providing good 
parameter estimates for forecasting. Since phytoplankton are largely affected by spatial inhomo- 
geneities in the ocean such as nitrogen runoff, regions of hypoxia, or upwelling, to name a few 
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parameter inhomogeneity- inducing effects, we wish to allow model parameters to vary spatially. 
These considerations are especially important since our models will be built over coastal domains 
where large changes in ocean biology occur spatially, leading naturally to spatially dynamic param- 
eters. 

We are only able to observe time series data Eq (4) as a movie and we know the model form 
of Eq (7), but want to estimate the parameters used to create the observed data. The system of 
Eq (7) will be taken as the drive system and we form a response system to be synchronized to the 
observations as, 



BP P7 

°— = AP + P(l-P)- 7 - h + K (P-P), (8) 

^ = AZ + kJ-^- -mZ + k(Z - Z), 
dt p + h V ' 

where we assume P(x, y, 0) ^ P(x, y, 0), Z(x, y, 0) 7^ Z{x, y, 0), k(x, y, 0) / k(x, y), and rh(x, y, 0) 7^ 
m(x,y). Thus, we do not know the initial model states, and wish to recover the spatially varying 
parameters m(x,y) and k(x,y). To derive Eq (8), a diffusive coupling term is added to each 
equation in Eq (7) accounting for the error between the drive and response values with a coupling 
strength, k. These additional terms drive P — > P and Z — > Z, so that the PDEs will synchronize 
after a short time. The synchronization is of identical type and dependent upon the choice of k, as 
is the synchronization speed. 



3 Results and Simulations of Autosynchronization Parameter Es- 
timation 

We modify the system Eq (7) as found in [5] by allowing the parameters to be nonnegative C°(£l) 
functions. Here O is the domain, which in the case of our simulations, C R 2 is a compact domain 
such as a rectangle or even a domain shaped as the Gulf of Mexico. Parameters are updated as 
diffusively coupled PDEs during the synchronization process as 

s(P -P) s > 0, (9) 
s(Z - Z), 

where we choose s = 30 for specificity and for which we observe good convergence results. The 
parameter equations are evolved simultaneously with Eq (8) with a forward Euler discretization 
and the same time step. The model form of Eq (9) was chosen after testing several forms and there 
may exist other forms for which synchronization is possible. Once the model form was chosen, 
a parameter search was performed to find s = 30. As we vary s and k, autosynchronization 
can fail, a common situation with diffusively coupled systems. Parameters may be updated as 
reaction-diffusion PDEs, by adding a diffusion term, however we need to restrict parameters to be 
nonnegative C 2 (f2) functions and stability may be affected. To begin the simulation, parameters 
are initialized as the constant function 

k{x,y,0) = 10, (10) 
rh(x,y,0) = 10. (11) 



dk 

dt 
dm 

Ik 
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We evolve Eq (7) forward and count the model output as observed data. Initial conditions for the 
response system are P(x, y, 0) = and Z(x, y, 0) = 0. Furthermore, to avoid values outside the 
normal range of Eq(7), we enforce that 

fP : < P < 2 
P = < : P < , 
I 2 : P > 2 

during the simulation. 

First, we develop synthetic datasets with spatially varying parameters to challenge our methods. 
Spatially dependent parameters are chosen to be in the range given in [5] for spatially irregular 
behavior. Three different functional forms for the parameters are tested. First, we use a linearly 
varying parameter set 

h(x,y) = a(- + — + l) +b, (12) 

mi(x,y) = c(— + — + l)+d, 
\n m I 

where m and n represent the size of the domain shown in Figure ?? and a = 2,b = 0.14, c = 
0.6, d = 0.1, and I = 5. Thus, appropriate parameters are chosen to maintain m(x,y) and k(x,y) 
in the target range. Next, we define a Gaussian parameter function as 



(x-n/2) 2 (y-m/2) 2 
2 + 2^ 

\2 



k 2 (x,y) = ae V *" ^ 7, (13) 
7712(0;, y) = ce 



where a,b,m,n have the same values. For example, k2(x,y) is shown on the top half of (a) in 
Figure 2. Finally, we define 



ks(x,y) = a cos(bx + d) sin(6y) + s, (14) 
rnsixiV) = c cos(6x + d) sm(by) + t, 

where a = 0.2,6 = n/(m/2),c = 0.6, d = tt/2,s = 0.5, and t = 1.5, to test the quality of the 
autosynchronization scheme to resolve fine structures in model parameters. The surface produced 
by k%{x, y) is displayed on the top half of (b) in Figure 2. We observe solution data at every time step 
relative to the response system, Eq (8), and the parameter system, Eq (9), to drive (P, Z) —> (P, Z) 
and (m(x,y),k(x,y)) — > (m(x,y),k(x,y)) as t -> 00. For brevity, only the parameters defined by 
Eq (13) and Eq (14) are shown and compared with their estimated counterparts. 

We observe autosynchronization for each test set of parameters and the spatial inhomogeneities 
in each case are effectively resolved. In Figure (1), the globally averaged error between the drive 
and response PDEs (Phytoplankton density) has been driven to less than 1.0 x 10~ 10 and the 
globally averaged error between true and estimated parameters has been driven to below 1.0 x 10 -6 . 
Interestingly, the results show a change of convergence rate midway through the simulation. For 
these results, we choose k = .3625. 

In, Figure (2), the results of autosynchronization after 2000 seconds of simulation are shown 
where reconstructed parameters are compared with their true counterparts. Both plots demonstrate 
how effectively the parameters have been reconstructed. Similar results were found by testing 
parameters that vary spatially according to Eq (12). 
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Synchronization for Phytoplankton 
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Figure 1: Left: Globally averaged synchronization error between drive and response PDEs (Phy- 
toplankton) on a log scale. Right: Globally averaged synchronization error between drive and 
response parameters on a log scale. Both plots obtained using model parameters in Eq (13). 



4 Syncronization by Sampling Only One Species 

To this point, an important criticism of our work is that we need to sample both species to drive 
the response model and parameters. As mentioned above, our interest in autosynchronization for 
parameter estimation stems from work with ocean models for phytoplankton-zooplankton ecology. 
In fact, hyperspectral satellite imagery provides phytoplankton density inferences but provides no 
data for zooplankton. Certainly, parameter estimation using the response model above will fail. 
Even given correct model parameters, it is impossible to forecast the model since zooplankton 
initial conditions are not supplied. Our problem of interest requires that we somehow estimate 
zooplankton initial conditions based on phytoplankton observations. 

We find that, by a modification of Eq (8), it is possible to drive zooplankton density to its true 
state by sampling phytoplankton alone. This is a first demonstration of the possibility of simulating 
this system with only partial knowledge. As an added bonus we observe autosynchronization. 
Thus this technique gives us a tool to estimate parameters and to initialize a model for short term 
forecasts. The response model that drives these results is 



BP P7 

— = AP + P(l-P)-^— + K (P-P), (15) 
oi P + n 

= AZ + k mZ, 

dt P + h 

where the absence of hats denotes where observation data is coupled directly into the PDE. We 
use a combination of diffusive coupling and complete replacement coupling in the response PDE to 
observe autosynchronization. Note that zooplankton density is no longer observed in Eq (15). The 
parameter update equations are 

at = ,li - p - p 1 (16) 
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Figure 2: Autosynchronization of spatially dependent parameters at t = 2000. Drive (top) vs 
response (bottom). Left: Model parameters given by k2(x,y) (shown) and m2(x,y). Right: Model 
parameters given by kz(x,y) (shown) and 7713(2;, y). 



with s\ = 0.2, S2 = 0.6, and k = 0.6. Here, as above, the parameter equations are evolved 
simultaneously with the (15) using a forward Euler discretization. Figure (3) shows results obtained 
when autosynchronizing with Eq (15) at t = 4434, a substantially longer time epoch. Results 
obtained for the Gaussian form, Eq (13), for parameters show that both k(x,y,t) and m(x,y,t) 
converge, as before, to their true values k2(x,y) and 7712(0;, y). Here, all initial conditions are set to 
P(x, y, 0) = Z(x, y, 0) = k(x, y, 0) = m(x, y, 0) = 1. 



drive parameter k drive parameter M 




Figure 3: Autosynchronization of spatially dependent parameters at t = 4434. Pictured as drive 
(top) vs response (bottom) with £2(2;, y) vs k(x, y, 4434) on the left and 7772(2;, y) vs m(x, y, 4434) 
on the right. 



In Figure (4), globally averaged errors are shown to diminish over time as the synchronization 
scheme evolves. Both parameters converge to within about 1.0 x 10 -6 after 4434 seconds. Impor- 
tantly, we note zooplankton convergence to nearly 1.0 x 10~ 7 of ground truth. Therefore, we need 
not sample zooplankton to observe autosynchronization and we find the true zooplankton density 
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Autosynchronization for Gaussian Parameters Synchronization forZooplankton 




Time Time 

Figure 4: Left: Globally averaged synchronization error between drive and response parameters 
on a log scale. Right: Globally averaged synchronization error between drive and response PDEs 
(Zooplankton) on a log scale. Both results are obtained using Eq (15) with drive model parameters 
k2(x,y) and m,2{x,y). 



profile such that model simulations may be initialized. 

5 Conclusion 

In this paper, we have shown that it is possible to derive an autosynchronization scheme for a 
system of PDEs. We emphasize here the improvements we have made upon past synchronization 
methods in that we use autosynchronization as a means of parameter estimation of parameters that 
exist in a function space. We assume that we know the model form of the true observed system, 
but have no prior knowledge of the parameters. By sampling at every time step, we observed 
identical synchronization between the response and drive systems as described in [2]. As a first 
attempt, we have given a model form for adaptive parameters in the response system such that 
we observe identical synchronization between response model parameters and true parameters, or 
autosynchronization. Our techniques were implemented on a benchmark model and results converge 
to ground truth. Thus, autosynchronization is observed for PDEs with spatially homogeneous 
parameters. 

Next, we considered the same system of PDEs wherein the parameters were spatially dependent. 
We provided a scheme with which we observe autosychronization of spatially dependent parameters. 
We tested our results against several different functional forms for parameters and found the method 
to be robust. 

We markedly improved upon these results once more with an autosynchronization scheme that 
requires sampling of only one species (phytoplankton). We noted that in order to evolve a system 
of PDEs for forecasting, we need initial conditions for both species; this is a serious problem when 
dealing with remote sensing data with which we can only observe one of the species. This concern 
was addressed by providing a response system that autosynchronizes parameters and synchronizes 
zooplankton using only phytoplankton data. These methods are therefore plausible for use in 
remote sensing problems. 

As discussed above, synchronization schemes can be proven to work for a given range of coupling 
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parameters using, for example, Lyapunov functions. It remains to be shown why this scheme works 
on this system, and to perhaps derive autosynchronization model forms for a wider class of reaction- 
diffusion PDEs. 

A drawback of this technique with application to hyperspectral satellite data is that data may 
be noisy; this is where filtering techniques have a built-in advantage. Data may also be occluded 
because of cloud cover. Since there is no hope to synchronize PDEs without data, we require 
techniques to fill in that which is missing, for example, inpainting. 

Another problem with applying these techniques to satellite imagery is temporal data reso- 
lution. There may be several days between images and autosynchronization requires ample data 
observations. The need for frequent observables is perhaps the main drawback to this method. 
However, autosynchronization may be advantageous for parameter estimation or model building in 
situations where spatiotemporal data are abundant and especially where parameters are expected 
to vary spatially. 
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